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ABSTRACT 

Context. There is an active debate about whether the properties of comets as observed today are primordial or, alternatively, if they 
are a result of collisional evolution or other processes. 

Aims. We investigate the effects of collisions on a comet with a structure like 67P/Churyumov-Gerasimenko (67P). We develop scaling 
laws for the critical specific impact energies Reshape required for a significant shape alteration. These are then used in simulations of 
the combined dynamical and collisional evolution of comets in order to study the survival probability of a primordially formed object 
with a shape like 67P. Although the focus of this work is on a structure of this kind, the analysis is also performed for more generic 
bi-lobe shapes, for which we define the critical specific energy g bil . The simulation outcomes are also analyzed in terms of impact 
heating and the evolution of the porosity. 

Methods. The effects of impacts on comet 67P are studied using a state-of-the-art smooth particle hydrodynamics shock physics code. 
In the 3D simulations, a publicly available shape model of 67P is applied and a range of impact conditions and material properties 
are investigated. The resulting critical specific impact energy greshape (as well as Q bil f° r generic bi-lobe shapes) defines a minimal 
projectile size which is used to compute the number of shape-changing collisions in a set of dynamical simulations. These simulations 
follow the dispersion of the trans-Neptunian disk during the giant planet instability, the formation of a scattered disk, and produce 87 
objects that penetrate into the inner solar system with orbits consistent with the observed JFC population. The collisional evolution 
before the giant planet instability is not considered here. Hence, our study is conservative in its estimation of the number of collisions. 
Results. We find that in any scenario considered here, comet 67P would have experienced a significant number of shape-changing 
collisions, if it formed primordially. This is also the case for generic bi-lobe shapes. Our study also shows that impact heating is very 
localized and that collisionally processed bodies can still have a high porosity. 

Conclusions. Our study indicates that the observed bi-lobe structure of comet 67P may not be primordial, but might have originated 
in a rather recent event, possibly within the last 1 Gy. This may be the case for any kilometer-sized two-component cometary nuclei. 

Key words, comets: general - comets: individual: 67P/Churyumov-Gerasimenko - Kuiper belt: general - 
planets and satellites: formation 


1. Introduction 

Comets or their precursors formed in the outer planet region dur¬ 
ing the initial stages of solar system formation. They may have 
been assembled by hierarchical accretion (e.g. Weidenschilling 
1997; Windmark et al. 2012b,a; Kataoka et al. 2013) or, alterna¬ 
tively, were born big in gravitational instabilities (e.g. Youdin 
& Goodman 2005; Johansen et al. 2007; Cuzzi et al. 2010; 
Morbidelli et al. 2009), thereby bypassing the primary accre¬ 
tion phase entirely. Whether the cometary nuclei structures as 
observed today are pristine and preserve a record of their orig¬ 
inal accumulation, or are a result of later collisional or other 
processes is much debated (e.g. Weissman et al. 2004; Mumma 
et al. 1993; Sierks et al. 2015; Rickman et al. 2015; Morbidelli & 
Rickman 2015; Jutzi & Asphaug 2015; Davidsson et al. 2016). 
The shape, density, composition, and other properties of comet 
67P/Churyumov-Gerasimenko (67P) have been determined in 
detail by the European Space Agency’s Rosetta rendezvous 


mission (e.g. Sierks et al. 2015; Hassig et al. 2015; Capaccioni 
et al. 2015). Based on this data, it has been suggested that its 
structure is pristine and was formed in the early stages of the so¬ 
lar system (Massironi et al. 2015), possibly by low velocity ac¬ 
cretionary collisions (Jutzi & Asphaug 2015). What is less clear 
is whether or not a structure like comet 67P would have been 
able to survive until today. 

The collisional evolution of an object of the size of comet 
67P was studied by Morbidelli & Rickman (2015) using dy¬ 
namical models of the evolution of the early solar system. In 
the “standard model”, as defined by the so-called Nice model 
(Tsiganis et al. 2005), cometary nuclei, or their precursors, orig¬ 
inated from an initial trans-planetary disk of icy planetesimals 
with a lifetime of a few hundred Myr. In this concept, the trans¬ 
planetary disk formed in the infant stages of the solar system 
beyond the original orbits of all giant planets, which were ini¬ 
tially closer to the Sun. This disk may have given rise to both the 
Scattered Disk and the Oort cloud (Brasser & Morbidelli 2013) 
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and thus, it may once have been the repository for all the comets 
observed today. According to the standard assumption, the dis¬ 
persal of the disk coincided with the beginning of the so-called 
Late Heavy Bombardment (Gomes et al. 2005; Morbidelli et al. 
2012), and had a lifetime of about 450 Myr before it was dynam¬ 
ically dispersed. 

As shown in Morbidelli & Rickman (2015), it is clear that in 
this standard model, an object of the size of comet 67P would 
have experienced a high number of catastrophic collisions and 
thus could not have survived. However, it was also shown that 
in the (hypothetical) case that the dispersal of the disk occurred 
early, right after gas removal, the collisional evolution of km-size 
bodies ending in the Scattered Disk would have been less se¬ 
vere, and a fraction of these objects may have escaped all catas¬ 
trophic collisions. We also note that in alternative models (e.g. 
Davidsson et al. 2016), the total number of comets is consid¬ 
ered to be lower than previously thought. Therefore, the fate of 
cometary-sized objects appears to depend upon the details of the 
dynamical scenario considered. 

However, whether or not an object like comet 67P would 
have been able to survive until today does not only depend upon 
its dynamical evolution but even more so on the “strength” of the 
body. In other words, it is crucial to know the critical specific im¬ 
pact energy at which the shape and structure of such an object are 
altered significantly. Previous studies of the collisional evolution 
of comet 67P (Morbidelli & Rickman 2015) used scaling laws 
for catastrophic disruption energies that are based on idealized 
spherical, solid icy bodies (Benz & Asphaug 1999), which may 
not represent well the properties of a highly porous cometary nu¬ 
clei. It is well known that the impacts in highly porous material, 
given its dissipative properties, can lead to very different results 
compared to impacts involving solid materials (e.g. Housen & 
Holsapple 2003; Jutzi et al. 2008). Furthermore, complex shapes 
such as the one of 67P may already be substantially affected by 
relatively low energy, sub-catastrophic impacts. 

It has been suggested recently that rotational fission and re¬ 
configuration may be a dominant structural evolution process 
for short-period comet nuclei having a two-component structure 
with a volume ratio larger than ~0.2 (Hirabayashi et al. 2016). 
In this model, the fission-merging cycle would begin once a two- 
component body enters the inner solar system and significant 
changes in the rotation period occur. The final shape of the comet 
nuclei (e.g. 67P) as observed today would then be the result of 
the last merger in this cycle. In this context, it is important to also 
study the survival probability of more general two-component 
structures, as such structures are required for the fission-merging 
cycle to begin. 

In this paper, we consider both the dynamical evolution and 
the response to impacts of objects with a 67P-like shape as well 
as generic bi-lobe structures. This combined approach allows us 
to compute the expected number of shape-changing collisions 
for such objects, as well as the related survival probability and 
possible formation age. 

In the first part of the paper, we describe our modeling ap¬ 
proach to study the effects of impacts on comet 67P and generic 
bi-lobe shapes. In Sect. 2, we determine the specific energies 
Greshape required to change a 67P-like shape significantly, as well 
as the corresponding Gbii for reshaping generic bi-lobe objects. 
The catastrophic disruption threshold Q ^ for bodies of 67P size, 
with the same properties, is computed as well here. Using the 
result of this modeling, we develop scaling laws for Greshape, Gbii 
and Gd* Finally, the simulation outcomes are analyzed in terms 
of impact heating and the evolution of the porosity. 


In the second part of the paper, we first describe the details 
of the dynamical simulations used in this study and discuss the 
differences and the improvements with respect to the previous 
work by Morbidelli & Rickman (2015; Sect. 3). We then com¬ 
pute the average number of shape-changing collisions for a body 
with a 67P-like shape as well as for generic bi-lobe shapes, us¬ 
ing the corresponding scaling laws (Greshape and Gbii)- In Sect. 4, 
the uncertainties of our model as well as alternative models are 
discussed, followed by conclusions in Sect. 5. 

A scenario of the late formation of 67P-like (two-lobe) 
shapes by a new type of sub-catastrophic impacts is presented 
in a companion paper (Jutzi & Benz 2017, hereafter Paper II). 


2. The effects of impacts on bi-lobe structures 

Here, in a suite of 3D smooth particle hydrodynamics (SPH) 
code calculations, we compute the specific impact energy 
Greshape required to significantly change the shape of comet 67P 
as well as of generic bi-lobe structures. The catastrophic disrup¬ 
tion threshold Q ^ for spherical objects of the same mass is com¬ 
puted as well. We consider a range of material (strength) proper¬ 
ties and various impact conditions. The simulation outcomes are 
also analysed in terms of impact heating and the evolution of the 
porosity. 

2.1. Assumptions 

Cometary nuclei come apart easily due to tides (Asphaug 
& Benz 1994) and other gentle stresses (Boehnhardt 2004). 
Laboratory materials analysis (Skorov & Blum 2012), observa¬ 
tions of comet disruptions by tides (Asphaug & Benz 1994) or 
fragmentation through dynamic sublimation pressure (Steckloff 
et al. 2015), suggest a bulk strength of <10-100 Pa for these 
weakly consolidated bodies. On the other hand, a high compres¬ 
sive strength of surface layers on comet 67P (Biele et al. 2015) 
was found at 0.1-1 m scales. For our analysis of the overall sta¬ 
bility, this kind of small scale (<~10 m) strength is not relevant, 
as we are interested in the bulk properties. In our modeling, we 
thus consider bulk tensile strengths of up to 1 kPa. The corre¬ 
sponding values of cohesion and compressive strength are -an 
order of magnitude higher (see Sect. 2.2). 

The low bulk densities of comets indicate substantial poros¬ 
ity; in the case of comet 67P it is about 75% (e.g. Patzold et al. 
2016). In our modeling approach (Sect. 2.2) it is implicitly as¬ 
sumed that porosity is at small scales and the body is homoge¬ 
nous. In the case of comet 67P, recent gravity field analysis 
(Patzold et al. 2016) indicate that the interior of the nucleus is 
homogeneous (down to scales of -3 m) and constant in density 
on a global scale without large voids. This suggests our approach 
of modeling a homogenous interior is justified. 

2.2. Modeling approach 

The modeling approach used here has already been success¬ 
fully applied in a previous study to the regime of cometesimal 
collisions (Jutzi & Asphaug 2015). We use a parallel SPH im¬ 
pact code (Benz & Asphaug 1995; Nyffeler 2004; Jutzi et al. 
2008; Jutzi 2015) which includes self-gravity as well as material 
strength models. To avoid numerical rotational instabilities, the 
scheme suggested by Speith (2006) is used. 


A61, page 2 of 13 


M. Jutzi et al.: How primordial is the structure of comet 67P? 


Table 1 . Material parameters used in the simulations. 


Type 

P e (Pa) 

Ps (Pa) 

Pso (kg/m 3 ) 

Po (kg/m 3 ) 

Pcompact (kg/m ) 

a 0 

A (Pa) 

P 

To (Pa) 

Ft (Pa) 

Low strength 

10 2 

10 4 

910 

440 

1980 

4.5 

2.67 x 10 6 

0.55 

10 2 

10 1 

Medium str. (nominal) 

10 3 

10 5 

910 

440 

1980 

4.5 

2.67 x 10 6 

0.55 

10 3 

10 2 

High strength 

10 4 

10 6 

910 

440 

1980 

4.5 

2.67 x 10 6 

0.55 

10 4 

10 3 


Notes. Crush curve parameters P e and P s (Jutzi et al. 2008), density of matrix material p s0 , initial bulk density p 0 , density of the compacted 
material p CO mpact> initial distention or 0 , bulk modulus A, friction coefficient p, cohesion F 0 , average tensile strength F T . 
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porosity, the compressibility is limited not by the EOS but by 
the crush curve of the P-alpha model. The elastic wave speed 
c G for a porous aggregate body can be very low, of the order of 
c e - 0.1 kms -1 . To take this into account, we apply a reduced 
bulk modulus (leading term in the Tillotson EOS; see Table 1). 
The approach has the additional major advantage that the time- 
steps become large enough to carry out the simulations over 
many dynamical timescales. Different values of c e = 10-100 m/s 
are investigated. 

The relevant material parameters used in the simulations are 
indicated in Table 1 . 


2.3. Setup and initial conditions 

2.3.1. Impacts on comet 67P and generic bi-lobe shapes 


Fig. 1. Pressure-porosity relations (crush curve) used in the simulations 
for the three different sets of parameters (low, medium, high strength) 
as defined in Table 1. Also shown are the results from laboratory exper¬ 
iments dust agglomerates (Guttler et al. 2009) and ice pebbles (Lorek 
et al. 2016). 


In our modeling, we include an initial cohesion Fo > 0 and 
use a tensile fracture model (Benz & Asphaug 1995), using a 
range of parameters that lead to an average tensile strength Ft 
varying from -10 to -1000 Pa. We consider Ft = 100 Pa as the 
nominal case. To model fractured, granular material, a pressure 
dependent shear strength (friction) is included by using a stan¬ 
dard Drucker-Prager yield criterion (Jutzi 2015). As shown in 
Jutzi (2015) and Jutzi et al. (2015), granular flow problems are 
well reproduced using this method. 

Porosity is modeled using a P-alpha model (Jutzi et al. 2008) 
with a simple quadratic crush curve defined by the parameters 
P e , P s , po, Pso and ciq. We further introduce the density of the 
compacted material p co mpact = 1980 kg/m 3 , which is used de¬ 
fine the initial distention ct 0 = p CO mpact/Po = 4.5 corresponding 
to an initial porosity of 1 — 1 /or 0 ~ 78%, consistent with obser¬ 
vations (Sierks et al. 2015; Kofman et al. 2015; Patzold et al. 
2016). (We note that p s o in this model is a parameter determin¬ 
ing the form of the crush curve and does not correspond to the 
density of the fully compacted material). As an estimate of the 
compressive strength <x c = PJ2 is used. As shown in Fig. 1, 
the pressure-porosity relations resulting from these parameters 
(for low, medium and high strength; Table 1) covers very well 
the range of experimental curves for dust agglomerates (Guttler 
et al. 2009) and ice pebbles (Lorek et al. 2016). 

We apply a modified Tillotson equation of state (EOS; e.g. 
Melosh 1989) with parameters for water ice. It is adequate for 
modeling the collisions considered here, where the most impor¬ 
tant response is the solid compressibility. As long as there is 


To setup the target, we apply a publicly available shape model of 
comet 67P 1 , which defines the surface of the body. Three differ¬ 
ent sets of material parameters as indicated in Table 1 are used, 
corresponding to different target strength. 

To determine ^reshape for 67P-like shapes, we investigate a 
range of impact energies using a range of impactor sizes of 
R v = 100-300 m and varying the impact velocities. Target and 
impactor both have the same initial material properties; their ini¬ 
tial bulk density is po ~ 440 kg/m 3 . We only consider impacts 
into the smaller of the lobes of comet 67P. Two different impact 
geometries are investigated (Fig. 2). 

To determine the critical shape-changing impact energy gbit 
in the case of more general bi-lobe structures, we set up a target 
consisting of two overlapping ellipsoids (Fig. 6). Each ellipsoid 
has an axis ratio of 0.6. The volume ratio between the two com¬ 
ponents is -0.5 and the total mass of the body is M t = 1 x 10 13 kg. 
For these targets, we only use the nominal set of strength param¬ 
eters (Table 1) and an impactor size of R v = 100 m. 

The simulations are carried out using a moderately high res¬ 
olution of -3 x 10 5 SPH particles. 

2.3.2. Collisions among spherical bodies 

In addition to ^reshape and Qbii, we also investigate the critical 
specific energy for catastrophic disruption Q ^ of spherical bod¬ 
ies of the same mass and material properties as in the model of 
comet 67P. To do this we consider 3 different size ratios of pro¬ 
jectile and target (1:2; 1:4; 1:8) and varying impact velocities. 
The impact angle is fixed to 45°. 


1 http://sci.esa.int/rosetta/54728-shape-model-of- 
comet-67p/ 
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Geometry 1 
c e = 10 m/s 


Geometry 1 
c e = 100 m/s 


Geometry 2 
c e = 100 m/s 


Geometry 1 
c e = 100 m/s 
rotating 
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Fig. 2. Shape-changing collisions on comet 67P. We use SPH to simulate impacts of a R p = 100 m projectile on the smaller of the two lobes of 
comet 67P. The minimal specific energy required to cause a significant change of the comet’s shape by such impacts, <2resha P e> is estimated for 
different impact geometries and rotation axis. The material strength is the same in all cases shown here (F T = 10 Pa). The effect of the material’s 
sound speed is investigated as well (top row ; in this case, a bulk modulus of A = 2.67 x 10 4 Pa instead of the nominal A = 2.67 x 10 6 Pa is used). 
Plotted is a surface of constant density which represents the surface of the comet; shown in red are regions on the surface with materials whose 
prescribed tensile strength was exceeded. As a rough average, the minimal specific energy required to cause a significant shape change is estimated 
as ^reshape ~ 0.2 + 0.1 J/kg, as indicated by the horizontal yellow line. 


2.4. Results 

2.4.1. Critical specific energy for shape change 

The results of our modeling of impacts on 67P are displayed 
in Figs. 2-5. We find that this particular structure, with two 
lobes connected by a neck, is significantly altered even by rela¬ 
tively low energy impacts. For a fixed set of material parameters 
(i.e. constant strength), the different impact geometries and ro¬ 
tation states considered here lead to slightly different outcomes 
(Fig. 2), but there are no major, order of magnitude, differences 
between the various runs. 

As it can be observed in Fig. 3, higher material strength lead 
to higher specific impact energy required to reach a certain de¬ 
gree of change in the overall shape. 

There is no unique way to define the critical shape-changing 
specific impact energy from these results, but rough estimates 
are possible. Based on visual inspection, we define Greshape for 
the different strength as: Greshape ~ 0.2 ±0.1 J/kg for F T = 10 Pa; 
Greshape ~ 1.0 + 0.5 J/kg for Yj = 100 Pa; Greshape ~ 2.0 ± 
1.0 J/kg for F t = 1000 Pa (Fig. 3) for the impacts with the R p = 
100 m projectile. For the simulations with the larger projectiles 
we obtain Greshape ~ 0.3 + 0.15 J/kg (R v = 200 m; Fig. 4) and 
Greshape ~ 0.15 + 0.075 J/kg (R p = 300 m; Fig. 5), using the 
nominal strength of F T = 100 Pa. These values are plotted in 
Fig. 7 and compared to the catastrophic disruption threshold, as 
discussed below. We note that impacts into the larger lobe may 


lead to slightly different values for Greshape? but we do not expect 
order of magnitude differences. 

The results of our modeling of impacts on generic bi-lobe 
shapes (using nominal strength properties) are displayed in 
Fig. 6. The estimated minimal specific impact energies for re¬ 
shaping are Gbii ~ 2.0 +1.0 J/kg, which is slightly higher than 
in the case of the 67P-like shape with the same strength (Gbii 
[F t = 100 Pa] ~ corresponds to Greshape for the F T = 1000 Pa 
case). 

2.4.2. Catastrophic disruption threshold 

The results of our modeling of catastrophic disruptions of spher¬ 
ical bodies with the same mass and material properties as in the 
model of comet 67P are shown in Fig. 7. We define the specific 
impact energy as 

Q = 0.5n r V 2 /(M t + M p ) (1) 

where /i v = M p M t /(M t + M p ) is the reduced mass, M p is the mass 
of the projectile and V the impact velocity. For the oblique (45°) 
impacts considered here, we also take into account that only a 
part of the mass of the colliding bodies is interacting (Leinhardt 
& Stewart 2012), and compute the Gd values of the equivalent 
head-on collisions. 

As expected, the energy threshold for catastrophic disruption 
Gd » Greshape, by -two orders of magnitude. 
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Fig. 3. Same as Fig. 2 but for different material strength Y T of the target. c e = 100 m/s in all cases. The critical specific energies are: g res ha P e ~ 
0.2+ 0.1 J/kg for Y T = 10 Pa (corresponds to average in Fig. 2); g reshape ~ 1.0 ±0.5 J/kg for Y T = 100 Pa; g res hape ~ 2.0+1.0 J/kg for Y T = 1000 Pa. 



v (m/s) 

5 

10 

20 

30 

Q (J/Kg) 

0.02 

0.07 

0.29 

0.66 


Fig. 4. Same as Fig. 3 but for R p = 200 m (F T = 100 Pa), g res hape ~ 
0.3 + 0.15 J/kg. 
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Q (J/Kg) 

0.06 
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1.00 

2.24 


Fig. 5. Same as Fig. 3 but for R p = 300 m (F T = 100 Pa), g res hape ~ 
0.15 + 0.075 J/kg. 
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Fig. 6. Results of impacts on generic bi-lobe shapes with nominal 
strength properties (Y T = 100 Pa) for two different impact geometries. 
R p = 100 m. The minimal specific energy required to cause a significant 
shape change is estimated as g b ii ~ 2.0 + 1.0 J/kg. 


Fig. 7. Critical specific impact energies gent- The energy thresholds for 
shape-changing impacts on a 67P-like shape (greshape) for different ma¬ 
terial strength are shown, as well as the catastrophic disruption energies 
<2p for various impact velocities. We note that the g b n values found 
for shape-changing collisions on generic bi-lobe shapes overlap the re¬ 
sults for greshape with Y T = 1000 Pa; they are not shown seperately. 
The solid lines show the scaling law (Eq. (2)) with parameters given 
in Table 2. The maximal global temperature increase d T shown on the 
right y axis is estimated by assuming that all kinetic impact energy is 
converted into internal energy: d T = g crit /c p where a constant heat ca¬ 
pacity c p = 100 J/kg/K is used. 


As found in previous studies (e.g. Jutzi 2015), in the disrup¬ 
tion regime, the results for g^ are almost independent of the 
material (tensile) strength. 

Our values of g^ for different impact velocities (Fig. 7) 
agree well with scaling law predictions (Housen & Holsapple 
1990), adopting a value for the coupling parameter of fi - 0.42, 
which is typical for porous materials. 


The g^ values for the weak, highly porous bodies considered 
here are slightly higher than the specific energies g^ ice found for 
solid bodies made of strong ice Benz & Asphaug (1999; Fig. 8). 
This result reflects the dissipative properties of material porosity 
and is consistent with previous studies (e.g. Jutzi et al. 2010). 

Also shown in Fig. 8 is the value of g^ suggested by 
Leinhardt & Stewart (2009) for weak ice as well as g^ predicted 
from scaling laws for collisions between gravity-dominated 
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IK 


Q for bodies of similar size colliding at 40 m/s (Davidsson et al., 2016) 
Q for bodies of similar size colliding at 1 m/s (Jutzi&Asphaug, 2015) 

Q D (Benz&Asphaug, 1999; solid ice, p 0 = 440 kg/m 3 ) 

Q D (Leinhardt&Stewart, 2009; weak ice, p 0 = 440 kg/m 3 ) 

Q* d (Leinhardt&Stewart 2012; "small bodies", c* = 5, p = 0.37) 

Q* d (this study; p = 0.42) 

Qreshape ( this stud y: P = °- 42 > Y T = 10 Pa) 

Qreshape ( this stud y: P = °- 42 > Y T = 100 Pa) 

Qreshape (this study; p = 0.42, Y T = 1000 Pa) 


10 



v (m/s) 


Fig. 8. Comparison of critical specific impact energies g crit . The scaling 
laws shown Fig. 7 are compared here with Q ^ values found in previous 
studies (Benz & Asphaug 1999; Leinhardt & Stewart 2009, 2012). Also 
displayed are the specific energies Q of collisions involving bodies of 
similar size (mass ratio of 1:2) for the bi-lobe forming collisions in study 
by Jutzi & Asphaug (2015) with very low velocities (V - Vesc ~ 1 m/s) 
as well as for collisions with a velocity of V = 40 m/s, corresponding 
to the average random velocity in the primordial disk during the first 
25 Myr in the model by Davidsson et al. (2016). 


bodies (Leinhardt & Stewart 2012). In these studies, the effects 
of material porosity were not taken into account. 

Finally, we also display in Fig. 8 the specific energies Q in¬ 
volved in collisions of bodies of similar size (mass ratio 1:2) in 
the bi-lobe forming low-velocity regime investigated by Jutzi & 
Asphaug (2015). As expected, those low-velocity (V - Vesc) ac¬ 
cretionary collisions have specific impact energies far below the 
disruption threshold. For reference, we also show the specific 
energy for collisions with much higher velocities (v = 40 m/s), 
which correspond to the average random velocity in the initial 
primordial disk in the model by Davidsson et al. (2016). For a 
mass ratio of 1:2, the specific impact energies are even above 
energy threshold for catastrophic disruptions Q 

2.5. Scaling laws for critical specific energies 

The results obtained in the previous section allow us to derive 
a go scaling law for porous cometary nuclei, which is a func¬ 
tion of impact velocity V and target size R (Housen & Holsapple 
1990): 

Q d = aR^V 2 -^ (2) 

where p and a are scaling parameters. 

For ^reshape and Qbii» we use a fixed target size R = 1800 m. 
As shown in Fig. 7, p = 0.42 also reproduces well the velocity 
dependence of these critical specific energies. The scaling pa¬ 
rameters for 2^, 2 re shape and 2bii are given in Table 2. 

2.6. Impact heating 

The effects of the impacts considered in this study (shape¬ 
changing impacts as well as catastrophic disruptions) are ana¬ 
lyzed in terms of impact heating and porosity evolution (below). 


Table 2. Parameters (SI units) for the scaling law 2crit = aR^V 2-3^, 
where R is the target radius and V the impact velocity. 


Scaling 

P 

a 

Qb 

0.42 

4.0e-4 

(2reshape (19 Pel) 

0.42 

9.0e-7 

Qreshape (100 Pa; nominal) 

0.42 

2.5e-6 

Qreshape (1000 Pa) 

0.42 

3.8e-6 

Qbii (nominal) 

0.42 

3.8e-6 


Notes. The scaling for shape-changing impacts on 67P (2resha P e) and 
for impacts on generic bi-lobe shapes (Q b u) on ly hold for a fixed size 
(R= 1800 m). 



Fig. 9. Cumulative post-impact temperature increase d T for two specific 
cases of shape-changing collisions, as indicated in the plot. 


First, in order to get an idea of the maximal global heating, we 
simply convert the total specific impact energy into a global tem¬ 
perature increase d T = 2crit/Cp where a constant heat capac¬ 
ity c p = 100 J/kg/K is used. The value of c p is a rough mass 
weighted average of the heat capacity of ice (Klinger 1981) 
and silicates (Robie & Hemingway 1982) at low temperatures 
(-30 K), assuming a dust-to-ice mass ratio of 4 (Rotundi et al. 
2015). Figures 7 and 8 (scale on the right) shows these d T values 
corresponding to collisions with a given specific impact energy. 

From this simple estimation, it is already obvious that im¬ 
pacts with energies comparable to 2reshape> the maximal global 
temperature increase must be limited to small values (d T <^c 
1 K). On the other hand, catastrophic impacts at kilometer scales 
may have the potential to lead to significant large scale heating, 
depending on how the impact energy is distributed. 

We compute the actual post-impact d T distribution for a few 
specific cases of the shape-changing (Sect. 2.4.1) as well as 
catastrophic collisions (Sect. 2.4.2). In the later, we only con¬ 
sider the material which ends up in the largest remnant (-50% 
of the initial target mass). The cumulative temperature distribu¬ 
tion in the case of the shape-changing impacts (Fig. 9) confirms 
that only a very small fraction of the material experiences sig¬ 
nificant heating. For the catastrophic collisions we find that the 
part of the target which experiences the largest impact effects is 
ejected. As a result, the material which is bound to the largest 
remnant (consisting -50% of the original target mass) is not af¬ 
fected much be the collision (Fig. 10) and the heating is limited 
(<1% of the mass is heated by dT > 1 K), even at relatively high 
collision velocities (600 m/s). 
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Fig. 10. Fraction of material in the final body that experienced a tem¬ 
perature increase larger than a certain value d T in catastrophic dis¬ 
ruptions with different velocities V. The mass of the largest remnant 
M\r/M t ~ 50%. Only the material belonging (i.e. which is bound) to the 
largest remnant is considered in the analysis. 

2.7. Porosity evolution 

Porosity is changed by impacts in multiple ways. First, material 
is compacted due to the pressure wave generated by the impact. 
On the other hand, material is ejected and the process of reac¬ 
cumulation of the gravitationally bound material can give rise to 
additional macroporosity. Our porosity model computes the de¬ 
gree of compaction (change of the distention variable). In order 
to specify the increase of the macroporosity, we treat each SPH 
particle individually according to its ejection/reaccumulation 
history. Particles which are lifted off the surface or are ejected 
and reaccumulated experience a density decrease, resulting in an 
increase of porosity. We assume that reaccumulated material can 
lead to the addition of macroporosity of maximal 40%, a typical 
porosity of rubble-pile asteroids (Fujiwara et al. 2006). To com¬ 
pute the total final porosity 0 to tai resulting from compaction and 
reaccumulation, we use the relation 

0total = 1 — 1/t^total (3) 

and define the distention 

total = min(pcompact/Pmiri5 ^max ) (4) 

where p m i n is the minimal density reached by the SPH particle 
and Pcompact = 1980 kg/m 3 . For this calculation we consider 
all particles which are gravitationally bound to the main body 
(largest remnant). The upper limit of the distention is given by 

t^max = <^0 CT V (5) 

where a v is the distention value corresponding to 40% macrop¬ 
orosity, a v = (1 - (p v )~ l with = 0.4, and ao = 4.5 is the initial 
distention. 

The resulting cumulative porosity distributions are calcu¬ 
lated for the same cases of shape-changing and catastrophic col¬ 
lisions as discussed in the previous section (Figs. 11 and 12). In 
the case of the shape-changing collisions, compaction is quite 
limited, even though the impacted lobe is severely disrupted. 
Because of the very low gravity, material is lifted off the sur¬ 
face/ejected by the impact. Due to the addition of macroporos¬ 
ity resulting from reaccumulation, the final average porosity is 
about the same as the initial porosity (Fig. 11). 





f 





Y T = 100 Pa, FL = 100 m, v = 200 m/s, includir 
" Y T =100Pa, R = 100 m, v = 200 m/s 

Y t = 100 Pa, FL = 300 m, v = 20 m/s, includir 
Y T = 100 Pa, R p = 300 m, v = 20 m/s 

ig macro-porosity - : 

, compaction only . " 

ig macro-porosity . 

, compaction only . 

Initial porosity - | 


0.001 0.01 0.1 1 
Mass fraction 

Fig. 11. Post-impact porosity distribution for two specific cases of 
shape-changing collisions, as indicated in the plot. The porosity 
calculation takes into account compaction as well as the increase 
of macroporosity. For comparison, the porosity distributions resulting 
from compaction only are shown as well. The final average poros¬ 
ity (compaction plus addition of macroporosity by reaccumulation) is 
78.8% (R p = 100 m) and 79.2% (R p = 300 m), respectively, while the 
initial porosity was 77.8%. 
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Fig. 12. Same as Fig. 11 but for catastrophic collisions. M\ r /M tot ~ 50%; 
only the material bound to the largest fragment M ]r is considered. The 
final average porosity (compaction plus addition of macroporosity by 
reaccumulation) is 83.3% (V = 10 m/s), 83.3% (V = 100 m/s), 83.3% 
(V = 200 m/s), 82.4% (V = 800 m/s), respectively, while the initial 
porosity was 77.8%. 


In the catastrophic disruptions, most of the material which 
undergoes collisional induced compaction does not remain on 
the main body (largest remnant). As a result, only ~10% of the 
material in the final main body has experienced significant com¬ 
paction. On the other hand reaccumulation plays a major role in 
this collision regime, resulting in a significant increase of macro¬ 
porosity. The final porosity is therefore even slightly higher than 
the initial porosity (Fig. 12). 

In Paper II, the interior porosity distribution of bi-lobe struc¬ 
tures resulting from sub-catastrophic collisions are compared to 
observations of comet 67P. 
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3. The combined dynamical and collisional 
evolution of comet 67P 

3.1. Modeling approach 

We follow the approach described in Morbidelli & Rickman 
(2015) in order to combine the dynamical evolution of the plan- 
etesimals precursors of Jupiter family comets with their colli¬ 
sional evolution. We do not repeat here a detailed description of 
the procedure, but we report on the differences and the improve¬ 
ments in the new implementation. 

These are of three kinds. First, we consider here only the 
dynamical dispersal of the original trans-Neptunian disk of plan- 
etesimals, which generates the Scattered Disk (the current source 
reservoir of JFCs). Thus, we neglect the phase ranging from the 
time when the gas was removed from the protoplanetary disk 
to the time when the giant planets developed a dynamical insta¬ 
bility that dispersed the planetesimal disk (Tsiganis et al. 2005; 
Gomes et al. 2005). This choice is made because Morbidelli & 
Rickman (2015) already showed that in the standard model, a 
comet the size of 67P has no chance to survive intact during this 
phase, if protracted for -400 My. On the other hand the debate 
on the timing of the giant planet instability is still open (see for 
instance Kaib & Chambers 2016; Toliou et al. 2016), so it might 
be possible that the aforementioned phase is short. There is no 
doubt, however, that the dispersal of the trans-Neptunian disk 
occurred and that this formed the Scattered Disk. In this case, 
Morbidelli & Rickman (2015) showed that during this process 
the number of catastrophic collisions for planetesimals the size 
of 67P is -1, so there might be some objects escaping break-up 
events. Thus, in this work we focus on this case, using improved 
assessments on the specific energies for catastrophic break-up 
and for reshaping, described in the previous sections. 

The second improvement over Morbidelli & Rickman (2015) 
concerns the dynamical simulations. Morbidelli & Rickman 
(2015) used the simulation of Gomes et al. (2005), which cov¬ 
ered only the first 350 My after the giant planet instability. This 
is when most of the action happens, but the subsequent 3.5- 
4.0 Gy cannot be neglected. Morbidelli & Rickman (2015) as¬ 
sumed that over this remaining time the orbital distribution of 
the Scattered Disk does not evolve any more, but its population 
decays exponentially down to 1 % of the original population after 
4 Gy. The 1% fraction comes from previous studies of the long 
term evolution of the Scattered Disk (Duncan & Levison 1997). 
Here we use the simulations presented in Brasser & Morbidelli 
(2013), which constitute a much more coherent set. Brasser & 
Morbidelli (2013) studied the dispersal of the trans-Neptunian 
planetesimal disk during the giant planet instability using a large 
number of particles (1 080000; including clones). At the end of 
the instability, they drove the giant planets towards their exact 
current orbits, so to avoid artefacts in the subsequent long-term 
evolution of the Scattered Disk. The evolution of the Scattered 
Disk was followed for 4 Gy. Because the number of active parti¬ 
cles decays over time, the test particles have been cloned 3 times, 
at 0.2, 1.0 and 3.5 Gy. In the final 0.5 Gy simulation, the parti¬ 
cles leaving the Scattered Disk to penetrate into the inner solar 
system as JFCs have been followed, in order to compare their 
orbital distribution with that of the observed comets. This final 
step is crucial to demonstrate that the Scattered Disk generated 
from the dispersal of the trans-Neptunian disk by the giant planet 
instability is a valid source of JFCs. 

The third improvement over Morbidelli & Rickman (2015) 
is that the collisional evolution is followed only for the parti¬ 
cles that eventually become JFCs in the final 0.5 Gy simulation. 



Fig. 13. Number of expected catastrophic collisions A disrup t during the 
formation and evolution of the Scattered Disk for the particles that even¬ 
tually become JFCs in the final 0.5 Gy simulation. Abrupt is computed 
using the scaling parameters for our new Q The symbols depict dif¬ 
ferent values for the exponent of the differential size distribution q , as 
labeled in the plot. 


These are 87 particles. We think that, potentially, this is an im¬ 
portant improvement. The particles that penetrate the inner solar 
system at the present time might have had specific orbital histo¬ 
ries relative to the other particles that either became JFCs much 
earlier or are still in the Scattered Disk today. Averaging the col¬ 
lisional histories of these three categories of particles may not 
be significant to address the specific case of 67P, which clearly 
became JFC only in recent time. 

Like in Morbidelli & Rickman (2015) the number of col¬ 
lisions suffered by each considered body is computed assum¬ 
ing that the initial disk particles represent a population of 2 x 
10 11 planetesimals with diameter D > 2.3 km, with a differential 
size distribution characterized by an exponent q. The minimum 
projectile size is determined by the scaling law (Eq. (2)) for the 
critical specific energy, with parameters given in Table 2. As for 
the exponent q, in agreement with Morbidelli & Rickman (2015) 
and previous studies of the comet size distribution, we consider 
here the cases with q = -2.5, -3.0 and -3.5. However, in the 
meantime the New Horizons mission to Pluto and Charon has 
measured the crater size frequency distribution, allowing the as¬ 
sessment of the size distribution of the trans-Neptunian objects 
larger than a few km in diameter (Singer et al. 2015). The pre¬ 
liminary results 2 suggest q = -3.3. Thus, we consider the re¬ 
sults for q = -3.0 and -3.5 as the most significant. However, we 
note that in alternative models (Davidsson et al. 2016) shallower 
slopes are preferred. 

We stress that the approach followed in Morbidelli & 
Rickman (2015) and in this work is conservative. This means 
that the number of collisions that are estimated is a lower bound 
of the real number of collisions. This is because the number of 
bodies assumed in the initial trans-Neptunian disk (2 x 10 11 with 
D > 2.3 km) is the minimum required, in absence of collisional 
comminution, to form an Oort cloud and a Scattered disk that 
contain enough objects to be sufficient sources of the LPC and 
JFC fluxes that we observe today. 


2 We note that based on the most recent results it has been suggested 
that there may be a deficit of small objects (Singer et al. 2016); see 
discussion in Sect. 4. 
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Fig. 14. Same as Fig. 13, but shown is the number of shape-changing 
collisions on a 67P-like body Af res hape> computed using scaling the pa¬ 
rameters for ^reshape for different strengths. We note that the number 
of shape-changing collisions A^n in the case of a generic bi-lobe shape 
with nominal strength properties is the same as 7V res hape for Ft = 1000 Pa 
(shown in the plot at the bottom). 


3.2. Results: number of disruptive and shape-changing 
collisions 

The number of events for each particle surviving in the Scattered 
Disk at the end of the disk dispersal simulation is shown in 
Figs. 13, 14. The results for the various types of collisions, using 
the corresponding scaling laws ( Q ^ and Greshape), are plotted. We 


Disruption, q: 
Disruption, q= 
Disruption, q: 
Shape-change, 1000 Pa, q= 
Shape-change, 1000 Pa, q= 
Shape-change, 1000 Pa, q= 


=-2.5 — 

Shape-change, 100 Pa, q=-2.5 


=-3.0 — 

Shape-change, 100 Pa, q=-3.0 - - ■ 


=-3.5 - 

Shape-change, 100 Pa, q=-3.5 


=-2.5 --- 

Shape-change, 10 Pa, q=-2.5 


=-3.0 — 

Shape-change, 10 Pa, q=-3.0 


=-3.5 - 

Shape-change, 10 Pa, q=-3.5 . 




Fig. 15. Cumulative fraction of particles (that eventually become JFCs) 
as a function of the number of collisions. This is an alternative represen¬ 
tation of the results already shown in Figs. 13 and 14. The solid lines 
correspond to the Q ^ scaling; the dotted lines were computed using 
Greshape (for different strength values). 



Fig. 16. Cumulative fraction of particles (that eventually become JFCs) 
as a function of the probability P(0) to escape all collisions. The dif¬ 
ferent line styles refer to different exponents for the differential size 
distribution q , as labeled on the plot. The three curves on the right cor¬ 
respond to the Gd scaling; the three curves on the left correspond to 
Greshape (with different strength values 10 Pa, 100 Pa and 1 kPa from 
left to right). For q = -3.0 and q = -3.5, the probability to miss all 
reshaping collisions is P(0) 10 -4 and the corresponding curves are 

not displayed here. 


note that the results for Gbii (impacts on generic bi-lobe shape 
using nominal material properties) are the same as in the case 
of Greshape with Yj = 1000 Pa (Table 2); they are therefore not 
displayed separately. 

Compared to the results by Morbidelli & Rickman (2015), 
the number of disruptive collisions is smaller (Fig. 13). This is 
mainly due to the fact the new Gd scaling law used here leads 
to disruption energies which are higher than the ones by Benz & 
Asphaug (1999) (which were used in the previous study). As dis¬ 
cussed in Sect. 2.4.2, this can be explained by the highly dissipa¬ 
tive properties of porous materials, which are taken into account 
in the new Gd- As Fig- 13 shows, for shallow size distributions, 
it is possible in principle that a significant fraction of the 67P 
sized objects escaped all catastrophic collisions. 
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Fig. 17. Mean number of reshaping collisions A res hape expected for 67P- 
like objects as a function of time for different strengths, as indicated 
on the y axis. We note that the number of shape-changing collisions 
Abu in the case of a generic bi-lobe shape with nominal strength prop¬ 
erties is the same as A^ res hape for Ft = 1000 Pa (bottom). Time t = 0 
corresponds to the beginning of the dynamical dispersal of the original 
trans-Neptunian disk of planetesimals, which generates the Scattered 
Disk; t - 4 x 10 9 yr is now. 


On the other hand, the number of shape-changing collisions 
(Fig. 14), requiring a much smaller impact energy (<2reshape), 
is substantially larger than the number of catastrophic events. 
As expected, the weaker the strength the larger the number 
of reshaping collisions taking place. Also, the steeper the size 




Fig. 18. Top: mean number of potential 67P-forming catastrophic colli¬ 
sions of a parent body of R = 3 km (computed using Q*^) as a function 
of time t (defined as in Fig. 17). Bottom : same, but for the scenario of 
67P formation by low energy sub-catastrophic collisions. 


distribution (larger q ), the larger the number of collisions hap¬ 
pening. However, in any case, even for the largest strength 
(1000 Pa) and the shallowest slope (q = -2.5), the number of 
reshaping collisions largely exceeds 1 for all comets. 

The results are summarized in Fig. 15 which shows the cu¬ 
mulative fraction of particles as a function of the number of col¬ 
lisions. In Fig. 16, the number of collisions N co \\ is converted into 
a probability to avoid all collisions F(0) = exp(-^V co n) and the 
normalized cumulative distribution of the P(0) values is plotted. 
The average number of collisions and the related probabilities 
are given in Table 3. 

It is also interesting to look at the number of collisions as a 
function of time (Figs. 17 and 18) as this in principle allows us to 
determine the time at which on average the last event of a certain 
type took place. 

For size distributions with q < -3.0, the last shape-changing 
event (on average) would have taken place in the last 1 Gy 
(Fig. 17), suggesting that the structure of comet 67P must have 
formed in a recent period. 

In Fig. 18, we plot the average number of events as a func¬ 
tion of time for two potential formation scenarios. In the first 
scenario, it is assumed that the structure of 67P formed as a re¬ 
sult of a catastrophic break-up of a parent body of R = 3 km. The 
corresponding number of collisions is then computed from our 
new <2^ scaling. In the second case, we consider impact energies 
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Table 3. Average number of shape-changing collision on a 67P-like object (A^ res hape). shape-changing collisions on a generic bi-lobe body (A^n) 
and catastrophic collisions (Abrupt). 


Type 

q = -2.5 

q = - 3 

q = -3.5 

N disrupt 

0.41 

(6.7E-1) 

0.79 

(4.5E-1) 

2.06 

(1.3E-1) 

^reshape (1C) Pa) 

4.92 

(7.3E-3) 

35.1 

(6.0E-16) 

258 

(7.3E-113) 

reshape (100 Pa) 

3.06 

(4.7E-2) 

18.1 

(1.4E-8) 

112 

(2.4E-49) 

N ]reshape (1000 Pa) 

2.53 

(8.0E-2) 

13.8 

(1.0E-6) 

79.6 

(2.7E-35) 

A^ii (nominal) 

2.53 

(8.0E-2) 

13.8 

(1.0E-6) 

79.6 

(2.7E-35) 


Notes. The corresponding probability P(0) to avoid all collisions is given in parenthesis. 


corresponding to a scenario of 67P formation by low energy sub- 
catastrophic collisions, as proposed in Paper II. Clearly, the num¬ 
ber of events in the later case are substantially larger. This sug¬ 
gests that it may be a more probable formation mechanism than 
the catastrophic break-up scenario (see a more detailed discus¬ 
sion on this topic in Paper II). 

4. Uncertainties and alternative models 

In this section we discuss some aspects of the robustness and 
uncertainties of our modeling approach and alternative models. 

4.1. Critical specific energies 

The values for the specific catastrophic disruption energies Q ^ 
are well defined and follow the expected scaling (Fig. 7). The 
critical specific impact energies for reshaping are not as well de¬ 
fined and do depend on the material properties. However, we 
explore a reasonably large range of material properties and also 
apply large error bars to the results in this case. In any case, there 
is no doubt that 2reshape ^ 2d and consequently, there must be 
many more shape-changing events than catastrophic disruptions. 

4.2. Dynamical model 

A crucial quantity in the dynamical model is the initial number 
of comets. The assumption of the existence of 2 x 10 11 comets 
is in line with estimates of the current Scattered Disk and Oort 
cloud populations and numerical estimates of the fractions of the 
primoridal disk that end up in these populations. Both could be 
wrong, in principle. However the fractions of the primordial disk 
population implanted in the Scattered Disk and Oort cloud that 
we use (from Brasser & Morbidelli 2013) are not very different 
from those found in quite different dynamical models (Dones 
et al. 2004 for the Oort cloud to Duncan; and Levison et al. 2008, 
for the Scattered Disk). Therefore, they seem to be robust. 

The number of comets used in our model are based also on a 
flux of Jupiter family comets which is assumed to be currently in 
a steady state. If this is not the case, the Scattered Disk could be 
less (or more) populated than predicted by the model. However, 
we find this unprobable for the following reason. The current es¬ 
timates for the populations in the Scattered Disk and the Oort 
cloud are consistent with these two reservoirs being generated 
from the same parent disk (Brasser & Morbidelli 2013). Thus, 
if the Jupiter family comet flux is now - say lOx - the mean 
flux (so to argue for a Scattered Disk lOx less populated), the 
same should apply for the flux of long period comets. But the 
fluxes out of Scattered Disk and Oort cloud follow different pro¬ 
cesses: for the Scattered Disk, this is resonant diffusion and scat¬ 
tering from Neptune; for the Oort cloud it is stellar perturbations. 


Therefore, it seems unlikely that both fluxes increased by the 
same amount relative to the mean values. 

Another crucial quantity in our modeling is the slope of the 
size distribution q , which determines the number of projectiles of 
a given size and thus the number of impacts with energies above 
the critical value. There is an ongoing debate about the form of 
the size distribution in the Scattered Disk population. We argue 
that the observations of the crater size distributions in the Pluto 
system by the New Horizons mission provides one of the best 
available constraints. The cratering of Pluto and Charon is dom¬ 
inated by the hot population (Greenstreet et al. 2015). All models 
agree that the hot population and the Scattered Disk population 
are the same population in terms of physical properties and ori¬ 
gin. In fact, the collisional evolution of the hot population is not 
more severe than that of the Scattered Disk. Both suffered most 
collisions during the dispersal of the primitive disk (or before, 
if the dispersal was late). It is true that comets have a shallower 
distribution (Snodgrass et al. 201 1) as well as have the craters on 
the Jovian satellites (Bierhaus 2006; Bierhaus et al. 2009). But 
this is probably because small comets disintegrate very quickly. 
On the satellites of Saturn, the crater size distribution is similar 
to the one expected from a projectile population with a size dis¬ 
tribution like that of the main asteroid belt (e.g. Plescia & Boyce 
1982; Neukum et al. 2005, 2006), i.e. it is the same as measured 
by New Horizons on Pluto and Charon. 

We note that based on the most recent analysis of the New 
Horizons data, it has been suggested (Singer et al. 2016) that 
the size distribution for small (<2 km) objects is shallower (q ^ 
-1.5) than at large scales. However, this result is still preliminary 
with uncertainties to be clarified. As discussed above, the TNO 
size distribution looks very similar to the size distribution of the 
asteroid belt, which is a result of a collisional equilibrium (below 
-100 km). This suggests that the size distribution of TNOs is in 
a collisional equilibrium as well. A change of slope below 2 km 
would produce waves in the TNO size distribution above 2 km. 
This is not observed, which may indicate that the change of slope 
is not as pronounced. 

To check the effects of a varying slope on our results, we per¬ 
formed additional calculations using q = - 3.3 for large (>2 km) 
and a shallower slope q s for small (<2 km) objects. We find 
that re-shaping collisions could be avoided for q s > -2, which 
means that if indeed q s = -1.5, a 67P-like shape would survive. 
However, we reiterate that this calculation considers a conserva¬ 
tive scenario without any collisional evolution in the primordial 
disk. 

4.3. Alternative models 

Alternative models to the standard model such as suggested 
by Davidsson et al. (2016) predict a much smaller collisional 
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evolution and are consistent with the idea of comets being 
primitive unprocessed objects, formed primordially. However, 
these models require the number of objects in the Scattered Disk 
to be orders of magnitude smaller. We note that there is no direct 
observational measure of the Scattered Disk population and all 
estimates are indirect and pass through models, so such a small 
number can in principal not be excluded. 

In is not clear, however, how bi-lobe structures would 
form/survive in these models. Previous studies indicate that the 
primordial formation of bi-lobed shapes, such as the one of 
comet 67P, by direct merging requires extremely low collision 
velocities of V/V esc ~ 1 (Jutzi & Asphaug 2015). This would 
have to take place at the very early stages of solar system forma¬ 
tion, probably while the gas was still present. In the later phases, 
relative velocities are much higher. In the model of Davidsson 
et al. (2016), average relative velocities are V = 40 m/s during 
the first 25 Myr. For kilometer-sized bodies this implies a ratio 
V/Vesc ~ 40! In fact, the corresponding specific impact energies 
are larger than the catastrophic disruption threshold (Fig. 8). Our 
results show that even relative velocities of a few m/s are destruc¬ 
tive and lead to reshaping (Figs. 3-5). Therefore, it is unlikely 
that primordial bi-lobe structures would survive this phase, and 
at the same time their formation by collisional merging is im¬ 
plausible due to the high relative velocities. 

5. Summary and conclusions 

We have estimated the number of shape-changing collisions for 
an object with a shape like comet 67P, considering a dynamical 
evolution path typical for a Jupiter family comet, using a “stan¬ 
dard model” of the early solar system dynamics. 

First, we computed the effects of impacts on comet 67P us¬ 
ing a state-of-the-art shock physics code, investigating range of 
impact conditions and material properties. We found that the 
shape of comet 67P, with two lobes connected by a neck, can 
be destroyed easily, even by impacts with a low specific impact 
energy. From these results, scaling laws for the specific energy 
required for a significant shape alteration (^reshape) were devel¬ 
oped. For more general applications, the critical specific energies 
to alter the shape of generic bi-lobe objects (Qb\\) was investi¬ 
gated as well. 

These scaling laws for ^reshape and gbii were then used to an¬ 
alyze the dynamical evolution of a 67P-like object and generic 
bi-lobe shapes in terms of shape-changing collisions. We con¬ 
sidered a conservative scenario without any collisional evolu¬ 
tion before the dynamical instability of the giant planets. Rather, 
we tracked the collisions during the dispersion of the trans- 
Neptunian disk caused by the giant planet instability, the for¬ 
mation of a scattered disk of objects and the penetration of tens 
of objects into the inner solar system. To do this we used a set 
of simulations (Brasser & Morbidelli 2013) that produces orbits 
consistent with the observed JFC population. 

We find that even in this conservative scenario, comet 67P 
would have experienced a significant number of shape-changing 
collisions, assuming that its structure formed primordially. For 
size distributions with q < -3.0, the last reshaping event (on 
average) would have taken place in the last 1 Gy. The prelimi¬ 
nary results of the New Horizons missions concerning the crater 
size-frequency distribution on Pluto and Charon suggest that the 
current trans-Neptunian population (i.e. including the Scattered 
Disk) has a differential power-law size distribution with an ex¬ 
ponent q ^ -3.3 (Singer et al. 2015). The possible consequences 
of a shallower slope for small (<2 km) objects, as suggested re¬ 
cently by Singer et al. (2016), are discussed in Sect. 4. 


It has recently been suggested that rotational fission and re¬ 
configuration may be a dominant structural evolution process for 
short-period comet nuclei with a two-component structure, pro¬ 
vided the volume ratio is larger than ~0.2 (Hirabayashi et al. 
2016). Our analysis of impacts on generic bi-lobe shapes shows 
that they would have experienced a substantial number of col¬ 
lisions with energies sufficient to destroy their two-component 
structure. This strongly suggests that the two-component body 
which is required to exist at the beginning of the fission-merging 
cycle cannot be primordial. 

Thus, according to our model, comets are not primordial in 
the sense that their current shape and structure did not form in the 
initial stages of the formation of the solar system. Rather, they 
evolve through the effect of collisions and the final shape is a 
result of the last major reshaping impact, possibly within the last 
1 Gy. A scenario of a late formation of 67P-like two-component 
structures is presented in Paper II. 

It is clear that the results presented here are based on the as¬ 
sumption that the standard model of dynamical evolution is cor¬ 
rect. Although some of its parameters are debated, as discussed 
in Sect. 4, we believe that the model is robust. We note that it 
is so far the only model which produces the correct number of 
objects in the inner solar system with orbits consistent with the 
observed JFC population. 

Our results clearly show that if this standard model of solar 
system dynamics is correct, it means that the cometary nuclei 
as they are observed today must be collisionally processed ob¬ 
jects. Therefore, the remaining important question is whether or 
not such collisionally processed bodies can still have primitive 
properties (i.e. high porosity, presence of supervolatiles). If this 
is not the case, then the standard model must be wrong. This 
would mean for instance that either the primordial disk was dy¬ 
namically cold and contained a much lower number of objects, 
as proposed by Davidsson et al. (2016) or that there is a lack of 
small comets, implying an abrupt change in the slope of the size 
distribution. 

However, the analysis of the outcomes of the detailed im¬ 
pact modeling carried out here (for shape-changing impacts 
and catastrophic disruptions) suggest that collisionally processed 
cometary nuclei can still have a high porosity, and could have 
retained their volatiles, since there is no significant large-scale 
heating. Therefore, they may still look primitive, meaning that 
the standard model is consistent with the observations of comet 
67R This question is investigated further in Paper II and also in 
an ongoing study of bi-lobe formation in large-scale catastrophic 
disruptions (Schwartz et al., in prep.). 

Primordial or not, the structure of comet 67P is an important 
probe of the dynamical history of small bodies. 
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